1. Field of the Invention
The present invention is directed to a 3D image reconstruction method, and in particular to such a reconstruction method suitable for use in computed tomography.
2. Description of the Prior Art
Computed tomography is a widespread diagnostic tool in medical technology and for material research. For example, in a two-dimensional version, a slice of the subject is transirradiated from different positions (directions), proceeding from an X-ray source. The attenuated intensities are then measured behind the subject by a one-dimensional detector array and electrical signals corresponding thereto are fed to a computer in digital form. The distribution of the attenuation values in the transirradiated subject slice is then determined by means of a mathematical reconstruction method in the computer. The most popular method is the filtered back-projection. In its general form, it can be characterized by the steps of (1) preprocessing (weighting) of the measurement values. (2) convolution of the measurement values (filtering), (3) weighting of the convoluted data, (4) back-projection, and (5) scaling of the result.
Above all, the convolution and the back-projection are computationally intensive and time-critical.
If the convolution is forgone, a very blurred image results--known as a layergram--which does not reveal any details.
In general, the discrete convolution is described by the following formula: ##EQU1## where x(n) represents the sequence of the measurement values, h(n) is the convolution kernel, and the result y(n) is the sequence of the convolved (filtered) data.
A value y(n) of the convolution result is thus a weighted sum of all the input values x(n). Since, in practice, a measurement value set x(n) only consists of a finite number, generally N, values, the above formula reduces to a finite sum. Thus, in order to calculate an individual convolved value, N multiplications and N additions are generally necessary. If, in the convolution kernel h(n) which is used, only M components are different from zero (i.e., have non-zero values), and if M&lt;N, then this leads to a further reduction of the computing outlay. Only M multiplications and additions, at the most, are then necessary for the calculation of a filtered value y(n).
The kernels utilized in computed tomography typically have full lengths; i.e., M=2N-1. The number of arithmetical operations is thus determined by N, the length of the measurement value fields x(n). These kernels are based on filter functions which are derived from ramp filters. The two most well-known are kernel based on the filter function formulated by Ramachandran and Lakshminarayanan, and the kernel based on the filter function formulated by Shepp and Logan.
Kernels utilized in CT typically have the following properties:
1. Symmetry, i.e. h(n)=h(-n). It is therefore sufficient to define only the central component h(0) and the right half h(n), for example, with n=1,2 . . .
2. Absence of zero-frequency component i.e., the sum over all kernel components is zero. This corresponds to the fact that in Fourier space, the ramp filter H(s) disappears at the zero point.
3. The central component h(0) is the single positive component; all others are negative or zero. Due to the final scaling of the total reconstruction result, the convolution kernel can normalized, without limiting its generality, such that the central component has the value h(0)=1. For n.fwdarw..infin., the absolute value of the kernel components goes toward zero.
The above described convolution in local space is equivalent to a component-by-component multiplication in Fourier space. After a certain kernel length, and given the utilization of kernels of full length, it can be numerically more effective to perform the filtering by multiplication in the Fourier space, i.e. to use the fast Fourier transformation (FFT) twice, and to use a very efficient multiplicative filtering in between: EQU FFT: x(n).fwdarw.X(s) EQU Y(s)=X(s)*H(s) EQU FFT (back): Y(s).fwdarw.y(n)
If the number M of the components h(n) which are different from zero is considerably lower than the length N of the measuring field, this leads to a reduction of the computing outlay only in the local space version, but not in the Fourier method.
A trend in medical technology is to survey and reconstruct an image of an entire volume of a subject in an optimally short time. This is achieved by the use of surface detectors instead of a one-dimensional detector array. A volume reconstruction from the measurement data obtained in a rotational angiography scan is also possible. FIG. 1 depicts the basic construction of such a known device. Such devices use the X-radiation significantly better than split-image tomography devices in which only a planar fan beam is emitted through a diaphragm. Such an angiography system (C-arm device) serves as exemplary embodiment for describing the background forming the basis of the invention.
The focus 1 of an X-ray source in the system shown in FIG. 1 travels around the subject 3 along an arc 2. An entire sub-volume of the subject 3 is transirradiated. The attenuated intensities behind the subject 3 are measured by a surface detector 4, usually an X-ray image intensifier. For this type of geometry, Feldkamp, et al., "Practical cone-beam algorithm;" J. Opt. Soc. Am., Vol 1, 1984, pages 1612-619, describes an efficient approximative 3-D reconstruction algorithm. This algorithm is also of the filtered back-projection type and strongly resembles the abovementioned 2-D algorithm. Specifically, in the Feldkamp algorithm the filtering of the measurement values of the surface detector is not performed in two-dimensional fashion, but only by rows, i.e. one-dimensionally, as in the classical split-image tomography with a single-row detector. All of the above discussion regarding convolution applies to this algorithm. From the pyramidal X-ray beam which emanates from the X-ray source focus 1, only the fan beam 5 is shown, and on the surface detector 4, only the detector row 6 is shown, formed by a series of detector elements.
Typical data for such a C-arm device consist of about 50 projections each with 1024.times.1024 measurement values. In the Feldkamp algorithm, there are over 5000 one-dimensional convolutions over measurement value fields x(n) of the length N=1024. Given the utilization of normal CT kernels of full length, such as that based on the Shepp and Logan filter function, this is an enormous computing outlay which leads to prohibitively long computing times for many applications, such as use during surgical or treatment intervention, for example.
In practice, in many cases the subject 3 is wider than the detector; i.e., the VOI (Volume of Interest) to be reconstructed is in the beam path, but the entire subject cross-section is not. FIG. 2 illustrates this situation. The measuring field is referenced 7 here. The measurement value profiles then no longer have the shape as in FIG. 3, with a drop-off to zero at the margins, but instead they exhibit sharp cuts, as depicted in FIG. 4. These unnatural sharp cuts are error points in the measurement data x(n), which propagate to all the data y(n) given the utilization of expanded kernels. Theoretically, an exact reconstruction is not possible from such split projections. In practice, however, truly usable results are achieved. Basically two methods are known to approach this problem, these being extrapolation of the measurement values, and transition to local convolution kernels.
An extreme example of a local kernel is a type known as the Laplace kernel, which has a double differentiation:
h(0)=1 PA1 h(1)=h(-1)=-1/2 PA1 h(n)=0 otherwise PA1 . . ,0,-1/4,-1/4,1,-1/4,-1/4,0, . . . or PA1 . . ,0,-1/6,-1/3,1,-1/3,-1/6,0, . . .
If this kernel is used, then the original subject f(x) is not reconstructed, but instead a modified distribution .lambda.f(x) is reconstructed. This has the same edges as f(x)--they are even exaggerated--but no longer gives information about the real value in homogenous regions. If the subject is a homogenous cylinder, for example, then the reconstructed values sag in the middle, as illustrated in FIG. 5a. Compared to the normal CT reconstruction, e.g. with the Shepp and Logan kernel, noise is transferred into the image center with higher intensity by the short Laplace kernel. To ameliorate this effect, smoothed versions of this differentiating core are also used, such as:
If the convolution is forgone, then the aforementioned layergram with a profile as in FIG. 5b arises. For mathematical reasons, the layergram is referenced .lambda..sup.-1 f(x). The omission of the convolution is equivalent to a convolution with a kernel known as the unit kernel, which is defined by h(0)=1 and h(n)=0 otherwise. The unit kernel does not have the typical characteristic of CT kernels that the sum over all components is zero. FIG. 5c depicts the exactly reconstructed profile of a homogenous cylinder.
A linear combination consisting of these two versions is designated as general .lambda.-reconstruction. The mathematical theory goes back to Smith et al., "Mathematical Foundations of Computed Tomography";Appl Optics, Vol. 24, No. 23, Dec. 1, 1985, pages 3950-3957. The basic idea is that, in the weighted superposition in originally homogenous subject regions, the concave and the convex behavior of the two individual versions at least partially compensate, and the margin information is nevertheless retained.
Another possibility is to cut off a standard CT kernel of full length, such as kernel based on the Shepp and Logan filter function, symmetrically with respect to the zero point after L&lt;N components, and to appropriately modify the remaining components such that the transition to zero occurs smoothly.
It should be noted that the goal of a reconstruction is not to achieve an optimally precise (as possible) computation of X-ray attenuation coefficients, but is to display diagnostically significant information for the physician. For this reason, different convolution kernels are used in medical computer tomography, dependent on the imaging task being addressed. The 3-D reconstruction from the data of a rotational angiography scan involves the displaying of fine, high-contrast vascular trees and vascular anomalies. FIG. 6 illustrates the boundary value field for a kernel optimization (region of the convolution optimization). The target parameters are, primarily, quality of result, and secondarily numerical efficiency.
Short kernels of the length M&lt;&lt;N have the properties of being local in their effect; i.e., possible points of error, as in spilt projections, remain local in their effect, and numerically efficient.